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The phase behavior of the primitive model of electrolytes is studied in the framework of var- 
ious mean field approximations obtained recently by means of methods pertaining to statistical 
field theory (CAILLOL, J.-M., 2004, J. Stat. Phys., 115, 1461). The role of the regularization of 
the Coulomb potential at short distances is discussed in details and the link with more traditional 
approximations of the theory of liquids is discussed. The values computed for the critical temper- 
atures, chemical potentials, and densities are compared with available Monte Carlo data and other 
theoretical predictions. 
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I. INTRODUCTION 



Various ionic systems including electrolyte solutions, molten salts, and colloids can be studied with a good approx- 
imation in the framework of the so-called primitive model (PM) which consists in a neutral mixture of M species 
of charged hard spheres (HS) of charges q a and diameters a a (a = 1, . . . , M) 1 . Of special interest is the restricted 
primitive model (RPM) where M = 2, the spheres have all the same diameter <r, and the cations and anions bear 
opposite charges ±q. In this paper we shall also consider the special primitive model (SPM) in which the number 
M of species as well as the charges q a are arbitrary but all the ions have the same diameter a. The case of a size 
assymetry of the ions will not be considered however. 

We have presented elsewhere an exact field theoretical representation of the PM obtained by performing a Kac- 
Siegert-Stratonovich-Hubbard-Edwards (KSSHE) transformation2iii2i£i§iL& of the Boltzmann's factor- ,10,11 . Thanks 
to this transformation, the grand-canonical partition function (GCPF) of the PM is essentially equal to the GCPF 
of a fluid of bare hard spheres in the presence of an imaginary random Gaussian held iip (tp real) with a covariance 
given by the Coulomb potential. Since the action H[ip] of the KSSHE field theory can be explicitely obtained then 
all the sophisticated techniques developed in statistical field theoryiSii 3 . can be applied a priori to the PM. Note 
ly-j ' that a KSSHE formalism can also be introduced to deal with the case of neutral fluidai^i^; the method is thus quite 
■^j- , general and has even been employed in numerical simulations^. Besides KSSHE formalism, various field-theoretic 
0^ ■ approaches have been proposed recently to study lattice and off-lattice versions of the PM. It is the place to cite 
the works by Di Caprio et alM^, the phenomenological (coarse-grained) field theory of Ciach and Stelli^aSSiSI, and 
■ recent applications of the method of collective variables by Patsahan and Mryglod22i2&2i. 

In ref. Ijlll) we have succeeded in obtaining the free energy of the homogeneous SPM at second order loop order in 



the framework of KSSHE field theory. This expression is used here to construct a Landau theory^ for the liquid- vapor 
(LV) transition of the model. Many mean-field (MF ) lik e theories of the LV transition of the RPM or the SPM have 
been proposed and are reviewed for instance in ref. I)2fj) . The aim of building MF theories similar to those discussed 
in this paper is not only to attempt to reproduce more or less accurately the LV coexistence curve of the RPM and its 
critical point, but also to provide a good starting point for a possible renormalization group analysis of ionic criticality. 
O ■ An interesting step in this direction was made recently in ref. l)23j) . 

Our paper is organized as follows. In sec. [H] we resume and extend the analysis of ref. lfTl|) on the KSSHE 
J> . representation of the SPM. Then we specialize to the case of a binary SPM in sec. 11111 where we give the expressions of 
the Landau function at the second order in the loop-expansion. This yields various MF theories for the LV transition 
which are studied numerically in sec. IIVI We conclude in sec. 



II. PROLEGOMENA 



A. The KSSHE transform 



We consider the three dimensional (3D) version of the SPM. The particles live in a domain Q C M 3 of volume |fi|. 
We work in the grand canonical (GC) ensemble. Let \x a be the chemical potential of the species a. For convenience 
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we shall define the dimensionless chemical potential v a to be v a — f3fi a (/3 = 1/kT, k Boltzmann's constant, T 
temperature) . 

Denoting by rf^ S the position of the i a -th ion of species a we note that only configurations lu = 
(JVi; rj, . . . , r]^ \ . . . \N a ;rf, . . . , \ ... \Nm', r^ 1 , . . . , rfJ M ) without overlaps of the spheres - i.e. such that \\r" — 
r ip\\ > cr - do contribute to the GCPF S[/3, {v a }, |^|]- It follows from this remark that the ions can be supposed 
to interact via a pair potential v a ^{r) = q a q0 w(r) where w(r) identifies of course with the Coulomb potential \/r 
outside the HS core, i.e. for r > er, but is otherwise arbitrary inside (i.e. for r < a). 

We have shown in previous works devoted either to neutral or charged fluids how to take advantage of the arbi- 
trariness of w(r) inside the HS core to define properly the KSSHE transforrr&ifiiiLii. In the case of the SPM the crux 
of the whole matter is to rewrite the electrostatic part of the configurational energy as a definite positive quadratic 
form. Therefore the regularization of w{r) for < r < a must be chosen in such a way that the Fourier transform 
w(k) is positive and w(0) is finite A2iiSiii Replacing the point charge q a of each ion by a radial distribution of charges 
q a T(r) smeared out inside its volume fulfills all these requirements and proves to be convenient. Thus, assuming 

r(r) = if r >a = a/2 , 

r(k) = 1 , (2.1) 
where r(k) denotes the Fourier transform of r(r) we indeed have 

= ^ > (2.2) 

as required and, of course w(r) = 1/r for r > a by virtue of Gauss theorem. Moreover for any reasonable distribution 
r(r) the self-energy q 2 a w(0)/2 of ion q a will be a finite quantity. For instance one can chose a uniform surface 
distribution of charges of radius a ~ a/2 <a, i.e. 

r a (r) = -^5(\\r\\-a), 

r a (k) = ^ , (2.3) 
ka 



yielding the regularized w a {r) 



w a (r) = — for r > a 
r 



2a — r 

= 5— for r < a . (2.4) 

a* 

Other types of regularization will be discussed in the remainder of the paper. 

Under these hypothesis a KSSHE transform can be performed and the GCPF of the SPM can be rewritten asii 

= (S HS [{v a + iq a 0}}) Va , (2.5) 

where z7 Q = v a + q 2 a w(0)/2 is a renormalized chemical potential and Shs [{^a + iqa4>}} is the GCPF of a mixture of 
bare hard spheres in the presence of the local chemical potentials v a + iq a (f>(f). The "smeared" field <p is defined to 
be the space convolution 



)=0V2 / dVarOln-raHMra) (2.6) 



which will be conveniently noted 



(1)=^t(1,%(2). (2.7) 



Note that in the whole paper, summation over repeated indices, either discrete or continuous, will always be meant 
(except if explicitly stated otherwise). 

The brackets (. . .) v in equation l|2.5|l denote a Gaussian average over the real scalar field ip(r), i.e. 



(...)„. = K. 1 fv<p...ew(-\(<p\v?\<p)) , 



(2.8) 
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where T>ip is the measure of the functional integration and 

W - 1 (l,2) = -i-A 1 5(l,2) (2.9) 

is the inverse of the positive operator v c (l, 2) = \jr\ 2 . We have also made use of the convenient notation 

(iplv^ 1 ^) = / d 3 ri / d 3 r 2 <p{r%) v^ 1 (r 1 ,r 2 ) <p(r 2 ) 
Jn Jn 

= ^(1) W - 1 (1,2)<^(2) . (2.10) 
After an integration by parts with appropriate boundary conditions one obtains the more transparent expressior&ii 

( v |« c -» = -L f d 3 r(V^) 2 . (2.11) 
To make some contact with statistical field theory we finally introduce the KSSHE action or effective Hamiltonian 

H[<P] = ~ (<p\Vc l \<P) - ^ S HS [{V a + iq a 4>}] , (2.12) 

which allows us to recast S under the form 

E[{u a }] = AC: 1 / V( f exp(-WM) . (2.13) 

It is important to distinguish carefully, besides the usual GC averages < ... >gcj between two types of statistical 
field averages : the already defined < . . . > Vc (cf equation H2.8[) ) and the < . . . >n defined as 

Since many thermodynamic quantities of interest can be expressed in terms of the charge correlation functions G^ 1 

(n) 

it is important to relate these functions to the field correlation functions G v . More precisely one has 

Gg°(l,...,n) = <pc(l)...po(n) >gc , (2.15a) 
GW(l,..,n) = < V (l)... V (n)> n , (2.15b) 

where, in equation H2.15a|) . the microscopic charge density pc is given by 

p C (l) = t{\,2) q a p a (2) , 

Mr) = ^<5 (3) (r-r.D, (2.16) 



the summation over the dummy indices a being understood in the first line of equation (|2.16() . It proves convenient 
to define also the truncated (or connected) n-body correlation functions by the relations 



•••>») = ^(1, ,n) - £ II G cV)^ — for n ^ 2 ' ( 2 - 1? ) 



where the sum of products is carried out over all possible partitions of the set (1, . . . ,n) into subsets of cardinality 

TH < tjSI. The relations between the truncated G^ T and G^ T were obtained in ref. fill) by a lengthy and cryptic 
method, a new simple derivation is given in appendix A. For a homogeneous system one finds at order n — the 
charge neutrality condition, i.e. 

PC = PaQa = , (2.18) 
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and, at higher orders, the following relations 

/?4 2) T (l, 2) = ^Ai*(l, 2) - ^JL_ Al A 2 Gf T (l, 2) , (2.19a) 
|3" /2 i"G«" )T (l,.,«) = i i ^A 1 ...A„G^ T (l,...,n) forrc>3. (2.19b) 

In ref. we have computed the free energy f({p a }) of the homogeneous SPM at the second loop order in the 
framework of the KSSHE field theory. Before giving and discussing this expression we want to stress that the exact 
f({Pa}) should be independent of the form of the pair potential w(r) inside the core, i.e. of the type of regularization 
adopted for the KSSHE transform. More generally, each term of say a systematic low fugacity or high temperature 
expansion of f({p a }) should also be independent of this regularization; this point was carefully checked in refs. l|9 HIo() 
in the case of the RPM. However here, as in ref. <fll|) . we consider a loop-wise expansion of f({p a })- The small 
parameter involved in this expansion cannot be given a clear physical interpretation and serves only to keep track of 
certain classes of Feynman diagrams^. Consequently each term of the loop expansion depends upon w(r) inside the 
core. Deciding which type of regularization to adopt is a matter of arbitrariness or mathematical skill. However it 
seems reasonable to impose that, at each order of the loop expansion, f({p a }) should be stationary with respect to 
the variations of w(r) inside the core. 

The zero-loop (mean-field (MF) or tree level), one-loop and two-loops expressions of f({p a }) will be denoted by 
f(°\ /W, and /( 2 ) respectively. Although the expressions obtained for f^, and f^ 2 ' in ref. ifTll) were derived 
with the assumption of a charge smearing regularization for w(r) one can safely assume that they remain valid for 
any reasonable regularization of w(r) inside the HS core (i.e. such that w(k) > and w(r) = 1/r for r > a). 



B. Zero-loop order 



For a homogeneous system one finds for the zero-loop free energy per unit volume the sloppy resulliii 



/3/ (0) (W) = flfas(«) " f Pall w(0) , (2.20) 

where fnsiiPa}) denotes the excess free energy per unit volume of the reference HS fluid. Moreover the charge 
neutrality condition p a q a = has to be be satisfied-^. Note that f3f^ diverges to — oo for a point like distribution 
r(r) = 5' 3 )(r) but remains finite after regularization if w(r) is well behaved a r — 0. Although disappointingly simple, 
equation (|2.20() can be however exploited for it can be shown that the MF free energy (3f^ ({p a }) constitutes an exact 
lower bound for the exact free energy /({p Q })±I. Maximizing the expression (HHt of f3f {0) ({p a }) with respect to he 
variations of w(r) inside the core should yield an optimized lower bound for f({p a })', unfortunately this mathematical 
problem has no solution since, as apparent from equation l|2.2U|l . /3f^({p a }) is a linear functional of w(r). However, 
if one restricts oneself to potentials of the form 12.21) the problem can be solvedii and one finds that the distribution 
r(r) which maximizes Pf^{{p a }) is a uniform surface distribution of charges of radius equal to that of the ions, i.e. 
the distribution r a (r) given by equation ^2. '61 with a = a , whence the optimized lower bound /3/hs({Po}) — l3pa<la/ a 
for the free energy, i.e. nothing but the well known Onsager boundiSS. We christen this regularization scheme of 
w(r) as the optimized mean field (OMF) scheme. 

The pair correlation functions h^^{r) = 9^^{r) — 1 and the direct correlation functions c^p(r)m at the zero-loop 
order are related to the free propagator A(r) of KSSHE field theory In Fourier space one has^i 

^S(fc) = hns, P {k) - t3q a q p A{k) , (2.21a) 
cL°i(*) = chs, P (k) - I3q a qpw(k) , (2.21b) 

where /ihs, p(k) and chs, p(k) denote the Fourier transforms of the ordinary and direct correlation functions of a HS 
fluid at the density p = J2 a Pa respectively. We emphasize that hf^r) and c^Ar), as given by equations (|2.21|) . do 
satisfy to the Ornstein-Zernicke (OZ) relations! 

Cm = c i%w + P^n * O) ( 2 - 22 ) 



5 



where the symbol denotes a convolution in space. The Fourier transform of the propagator has the following 
expression^ 

/Ufc) — w(k) ^ 

1 + P p a qi w{k) ' 

Therefore w(r) and A(r) also satisfy to an OZ-like equation in direct space, that is 

A(r)=w{r)-/3[p a q 2 a ]A*w(r). (2.24) 

We stress that the propagator A(r) as well as the correlation functions \(r) and cf'Jr) do depend on the form of 
w(r) inside the core. 

C. One-loop order 

The one-loop free energy per unit volume is given byii 

1 r ifllt 

f3f {1 \{ Pa }) = f3f (0 \{p a }) + 2 J (2^)3 Hl + P\P«&Mk)) 

1 f d 3 k 

= Pfm({ Pa }) + 2 J 7^3 t ln C 1 + P\P«&\*> (*)) ~ P[p<*ql]w (k)] ■ (2.25) 

/3/*- 1 - 1 obviously still depends on the pair potential w(r) inside the core but we note by passing that the expression (|2.25() 
of ({Pa}) remains finite (contrary to that of /3/*- ') for point-like distributions of charges, i.e. for tdhM = S^ 3 \r). 
One indeed obtains in this case the well-known Debye-Hiickel (DH) result 

/3/dh({p q }) = Pfim({Pa}) - , (2.26) 

where k 2 = Att/3 PaQa is the squared Debye number. Incidently the DH propagator can also be computed with the 
result 

A DH (r) =exp(- K r)/r . (2.27) 

In the Gaussian approximation the free energy is given by its one-loop expression l|2.25|) and the correlation functions 
by their zero-loop expressions (^2.21f! lx i 12 . The Gaussian approximation of the KSSHE field theory therefore coincides 
exactly with the usual random phase approximation (RPA) of the theory of liquids 1,29 . One can further demand 
the one-loop free energy /3/^({Pa}) to be independent of the expression of w(r) inside the core. As well known the 
stationary condition 

; ; — r = IOT T < O , 2.28) 

ow{r) 

yields the conditions g^p(r) for < r < a and for all pairs (a,/3). We recover in that way the optimized random 
phase approximation (ORPA) of the theory of liquid a 1 ! 29 . It follows from eqs. 12.21fl that the condition of nullity of 
5a°^( r ) ms ide the core can be rewritten as 

w(r) = 1/r for r > a 

A(r) = for r < a , (2.29) 

where w(r) and A(r) are linked by the integral equation l|2.24|l . Note that this equation does not involve the correlation 
functions of the reference HS system and that it also appears in the mean spherical approximation (MSA) of the SPM 
(in this case the exact correlation functions /ihs and chs in the right hand side of eqs. I|2.21|l have to be replaced by 
their Percus-Yevick (PY) expressions , i.e. hpy and cpy respectively^) . It turns out that equations l|2.24|l and (|2.29() 
can be solved analytically with the resulbi^ 

wmsaM = 1/r for r>cr, (2.30a) 
B Bt 

wmsa(0 = — (2 ) for < r < a , (2.30b) 

(7 a 

B = * 2 +*-*d + 2*) 1/2 f (2 .30 C ) 
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where x = kct. We remark that the function wmsa(0 meets all the requirements of a regularized potential, i.e. 
^msaM — l/ r f° r r > a and w(k) > as can be shown easily, and therefore can be used as a regularizator of the 
KSSHE transform. It must be stressed that, despite of the formal analogies between equations l|2.30bjl and (|2.4|l -the 
latter being deduced from the former by setting B = 1-, wmsa(^) cannot be interpreted as the interaction energy 
between two smeared distributions of charge since < B < 1 as discussed in ref. pij) . The same function u>msa(V) 
appears in the MSA and ORPA approximations for the SPM and yields in both cases ion-ion pair correlation functions 
which vanish inside the HS core. By contrast, the regularized potential w a {r) which enters the OMF free energy does 
not insure the nullity of g^pir) for r < a. 

With the choice w(r) = wmsa(t) the one-loop free energy (|2.25|l can be computed exactly. The result is well- 
known^i but we cannot refrain from giving the following derivation which seems to be original. Let us first compute 
the internal energy per unit volume worpa = d 'f3 /orpa({/9(i}) / 'd '(3 . It follows from equation l|2.25|l that 

1 f d 3 k ~ d 

WORPA = J J2^f 13 A msaO) [P [p a ql] WMSA(k)] 

= ~\j d 3 r(3[p a q 2 a ] AmsaW^ [P [p a ql] w MSA (r)] , (2.31) 

where we made use of the expression (|2.23ll of A(fc). Since the propagator Amsa('") vanishes inside the core in the 
ORPA (and MSA) approximations we have merely to compute the derivative of /? [Pct9„] wmsa(t) with respect to f3 
for r > er which is trivial since wmsaM = 1/r in this case. Therefore 

g r 

"orpa = -- [p a q 2 a ] I d3r Amsa(^) [p a q 2 a ] wmsa(^) 

= | [ Pa ql] (Amsa (r = 0) - wmsa (r = 0)) , (2.32) 
where we made use of equation l|2.24|l to obtain the last line. Since Amsa(0) = 0, we have finally 

—ft 

WORPA = — [p a q a ] WMSA(O) , 

= -p[p a ql]B/a , 

x 2 + x - x(l + 2x) 1/2 

47T/3CT 3 ' 

which leads, after integration, to the ORPA free energji22i2£ 

/3/orpa({/Oq}) = /3/hs({p q }) + dp mqrpa(/3 ) 



(2.33) 



(i 



at n ^ 3x2 + Gx + 2-2(1 + 2a;) 1 / 2 ,„ OA . 

Pfm({p a }) tx— t — ■ ( 2 - 34 ) 



127TCT 3 



D. Two-loop order 

The free energy of the SPM at the second-loop order has the following expressionil 

/3/ (2) (W) = Pf {1) ({p a })- ^lp a q 2 a } 2 Jd 3 rh HS , p (r)A 2 (r) 

+ ^WqlfJ d 3 rA 3 (r) , (2.35) 

where /jhs, p(f) denotes the usual pair correlation function of a fluid of a single species of hard spheres at the total 
number density p = ^ Q p a . Demanding the independence of 13 f 1 ^ ({p a }) with respect to the variations of w(r) inside 
the core is a formidable mathematical task that cannot be achieved analytically. Each of the regularization schemes 
that we have discussed previously can be considered however, yielding different expressions for the propagator A(r) 
and thus for (3f( 2 \ Expression l|2.35|l remains however tractable, at least numerically, since it involves integrals of the 
functions A(r) and /ins, p( r ) which can be both evaluated numerically (for instance, the approximation /ihs, p( r ) ~ 
^py, p{r) can be used safely at low and moderate densities). 



7 



E. The WCA regularization scheme 

We end this section by discussing briefly the popular Weeks-Chandler- Andersen (WCA) regularization scheme. In 
the WCA scheme the potential w(r) inside the core is assumed to be a constant equal to l/cs2i. Clearly it amounts 
to define w(r) as the electric potential created by a uniform surface distribution of charge T a {r) where T a (r) is the 
distribution defined at equation (|2.3|l with a = a. It follows from this remark that the Fourier transform of the 
regularized pair potential reads as 

i%ca(k) = — p — 

with r a = sin(ka) / (ka) . Therefore u5wca(&) is not positive for all k and the KSSHE transformation is ill-defined. As 
a consequence the propagator Awca(0 can be singular (cf equation l|2.23[l ). This singularity yields the so-called RPA 
catastrophe 1 - for the one- loop (RPA) free energy since the argument of the In in the RHS of equation l|2.25[l can become 
negative. This finite wave-number k instability, which however is absent in the OMF and MSA regularization schemes, 
has been made responsible of the possible existence of charge-ordered phasei2*22i2ii2i. Although an order-disorder 
transition of the lattice- version of the RPM was indeed observed in Monte Carlo simulations^SiS, the evidences of a 
similar transition for the off-lattice version of the model are lacking at the time of writing as far as the author is well 
informed. 



III. THE BINARY MIXTURE 



In this section we specialize the results of section [H] to the case of a binary mixture (M = 2). The cations and 
the anions have all the same diameter a and bear charges q\ = q and q2 = —£q respectively where £ is the charge 
assymetry parameter. It follows from the neutrality condition (|2.18|l that the ion densities are p\ = £p/(l + £) and 
P2 = p/(l + (,) where p = p\ + pi is the total number density. For convenience we introduce the reduced densities 
Pa = Pa& 3 and free energies per unit volume f3f = /3/er 3 . The reduced free energy /3/hs of the reference HS system 
is thus given by 

^7H S (^,P2) = /37Hs(P*)+^( T ^lne-ln(l + 0) , (3-1) 

where /3/ H (p*) is the free energy of a pure HS liquid at the reduced density p*. The additional term in the RHS of 
equation l|3.1|l corresponds to the mixing entropy contribution. Following the authors of ref. 1)341) we define the reduced 
inverse temperature f3* = (3q 2 ^/a and the reduced temperature as T* = 1/(3*. The reduced Debye number will of 
course be defined as n* = na (note that k* 2 — 47r/3*p*) and we will also frequently make use of the dimensionless 
squared Debye number A = k* 2 . 

Quite remarkably the electrostatic contribution to the one-loop free energy depends solely on parameter A and not 
on the charge asymmetry factor £. Indeed one has 

Jf {1) =A/ HS (P*)+P* (j-^ln^-lntl + o) +/[A] , (3.2) 

where the expression of the function I[X] depends on the choice of regularization made for w(r). The expressions of 
I[X] in the DH and ORPA schemes have been given at eqs. (|2.26|) and l|2.34|l respectively. In the OMF and WCA 
schemes one has 

2 f°° , 2 A A Asina; 4 \ Asin:r 4 \ 
W[A] = dxx 2 (h^l + I — —J - Z — -) , (3.3a) 

^wca[A] = ^— 2 J dxx 2 ^ln ^1 H — g-sina;^ ^sina;^ , (3.3b) 

respectively. Both integrals occurring in eqs (|3.3|l are convergent at x — and x — 00 and can be handled numerically 
easily. However, the second integral (|3.3bl) is singular for large values of A (the so-called RPA catastrophe) since the 
argument of the In can become negative at finite x. The stability limit is found to be A < A max . ~ 84.2. 

In the case of a regularization by a surface distribution of charges of radius a = a/2 < a one has simply / a [A] = 
I a [Xa 2 ]a~ 3 = IoMF[Aa 2 ]a~ 3 - It is easy to check that in the point-like limit a — > one gets I a [X] -> -A 3 / 2 /12, which 
coincides with the DH result (|2.26|) as expected. 
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At the two-loop order the free energy depends now not only on A but also on the density and the charge asymmetry 
parameter £. One deduces from equation l|2.35[) that 

Tf {2) =Jf (1) \K[\,p*] + ± { ±=^L[\,p*} , (3.4) 

with 

r°° _ 2 

K[X,p*] = / dx iirx 2 A x (x) hns, P (x) , (3.5a) 
Jo 

1 f°° 

L[X,p*} = — / dxAnx 2 Al(x) , (3.5b) 
P* Jo 

where x = r/a is the distance in reduced units and A^(x) = (3*p*a~ 2 A\(r — ax) is the reduced propagator. We 
stress once again that, since A\ is simply related to the pair potential w (cf equations l|2.23[) . and (|2.24[l ) it depends 
therefore on the regularization scheme. For a given regularization scheme it depends only upon the parameter A. 

(2) 

Moreover, we note that (3f is invariant under the transform £ — > l/£ as expected and that the second term in the 
RHS of equation 1)3. 4JI vanishes in the case of the RPM (£ = 1). Expressions 1)3. 5a|) and (|3.5b|) can easily be computed 
numerically for each regularization scheme (for numerical details see Appendix B). 



IV. THE LIQUID- VAPOR TRANSITION OF THE BINARY SPM 
A. Theoretical background 

In past recent years extensive Monte Carlo simulations have been performed to study the liquid-vapor transition 
of the RPM and to elucidate the nature of ionic criticality. It seems now well established that the critical exponents 
of the RPM are those of the 3D Ising model2&i2£i2LiS. The more recent estimates of the Orsay group for the 
critical temperature, density, and chemical potentials are^ T* = 0.04917 ± 0.00002, p* = 0.080 ± 0.005, and v* = 
— 13.600 ± 0.005 respectively (the critical pressure is not known with accurate precision). Recent numerical studies 
are also available for binary mixtures with an asymmetry in charge or/and in siz o 34 i 39 . 

Here we consider only the binary SPM and study various MF theories built with the approximate free energies 
discussed in section IIIII We work in the framework of the Landau theory with a Landau function given by the 
expression^ 

v)(y,P,<p)=pf(j3,<p)-v<p, (4.1) 

the minimum of which gives minus the MF pressure, i.e^ 

-#?WM) = waw(v,0,(p) . (4.2) 
v 

Note that ip plays the role of the order parameter of the Landau theory. The free energy function (3f{(3, tp) which 
enters equation (|4.1(l can be rewritten as 

(3f(P,<p)=0fM<P)-<p(j^to£-ln(l + ^+J(\,<p) , (4.3) 

J(A, <p) = \k{\, 9 ) + lii^L[A, ip] . (4.4) 

/3/hs(v 3 ) i s the free energy of the HS reference system as given by equation l|3.1|l and the functions /, K, and L in 
the RHS of equation (|4.4|l are those defined in section fTTTl Note that we have dropped all the "*" to simplify the 
notations, reduced quantities are implicitly meant. 

Above T c and for any chemical potential v , w(u, f3, ip) is a convex function of the order parameter ip and equation l|4.2|) 
admits a unique solution Tp~(y, (3) which is the (unique) solution of the stationary condition 



where 



dw 
dp 



= . (4.5) 





w 



q> 9 

FIG. 1: The Landau function w(p) and its first three derivatives w^(ip) at T < T c . The functions i// 2 ' and ii/ 3 ' are 
independent of the chemical potential. The zeros of m/ 2 ' give the spinodal densities, w^ 1 ' and w are diplayed for arbitrary v 
(dashed curves) and for v = v cocyi (solid curves). At the coexistence the grayish areas are equal. 



Therefore the MF pressure is given by (3Pmf{v, ft) — —w(v, ft, ip) which yields the expressions of the MF densities and 
energies per unit volume 

i q\ 9/3Pmf _ _ v 

Pmf{v,P) = — 5 =<p, (4.6a) 

av 

um?{v,P) = - 3 -^ E - = ^P » ( 4 -6b) 

where we made use of the stationary condition l|4.5(l . 

Below T c , w(y, ft, ip) is no longer a convex function of tp and a second minimum of w(v, ft, ip) arises for some range 
of chemical potentials. For T < T c (ft > ft c ) there is a unique chemical potential v CO c^{ft) corresponding to two values 
Tp g and Tp l > Tp g such that —ftPuF, coox(/3) = w(v cocx ,ft,p g ) = w(j/ CO cx, ft,Pi) (see figure (QJ). These two minima 
corresponds to the coexistence of a gas phase at the density pmf, g (ft) = tp g and a liquid at the density pmf, i (ft) = Pi ■ 
At the critical temperature the two solutions p g and p l merge into a unique solution Tp c = pmf. c corresponding to 
the critical density. 

We have determined numerically the coexistence curve of the SPM for the one-loop and two-loop expressions of the 
free energy for various regularizations schemes discussed in section [HJ i.e. 

• charge smearing regularization scheme (i.e. w(r) = w a (r)). We limited ourselves to the case of the surface 
distribution (|2.3H with a* = a/a ranging from (DH case) to 1 (OMF case). 

• MSA regularization scheme (i.e. w(r) = Wmsa(j*))- 



• WCA regularization scheme (i.e. w(r) = wwca(0j c f sec - C I n this case a numerical solution was obtained 
only for temperatures not to far below T c in order to avoid the RPA catastrophe. 




FIG. 2: Coexistence curves of the RPM in the one-loop approximation with w(r) — w a (r). From top to bottom a* = (DH), 
a* = 0.2, a* = 0.4, a* = 0.6, a* = 0.8, a* = 1 (OMF). Solid lines : coexistence curves, dashed line: spinodal curves. 

B. Numerical results 

We have found worthwhile to discuss shortly the algorithm devised to obtain the coexistence curve in appendix B. 
The results for the critical temperatures, chemical potentials, pressures, densities and energies are summarized in 
tables (fll ITT1 ITTTli and are discussed below. 

1, One-loop results 

We have already noticed that the free energy of the binary SPM is independent of the charge assymetry factor 
£ at the one-loop order. The coexistence curves and notably the critical parameters (T*, p*,.. .etc) are therefore 
independent of £ which is a serious drawback clearly in contradiction with numerical simulation resulta^Si. 

Our data for the one- loop approximation are reported in table (HJ. In the case of the MSA regularization scheme 
one recovers the MSA theory provided that the HS free energy is chosen to be that of the PY approximation. In 
this case we observed a perfect agreement of our data those obtained by Gonzales- Tova:»si. We stress that the results 
reported in table {IJ for the MSA regularization scheme were obtained, as all other results reported in all the tables 
of this paper, by making use of the Carnahan-Starling approximation for the HS free energy and of the PY solution 
hpy(r) of the HS corelation function. 

As apparent in table {j} and figures (J2J and (|3J) the coexistence curves depend strongly upon the regularization 
scheme adopted for the Coulomb potential inside the HS core. This is particularly striking for the charge smearing 
regularization scheme where the critical temperature T* decreases steadily from T* — 0.3271 for a* = (DH case) 
to T* = 0.1150 for a* = a/ a = 1 (OMF case). Similarly the critical density p* decreases with the diameter a of the 
charge distribution from p* c = 0.12267 for a* = to p* = 0.02198 for a* = 1. 

The coexistence curves at the one-loop order obtained with the OMF, WCA, and MSA (or ORPA) regularization 
schemes are displayed in figure The MSA scheme gives the better critical temperature, i.e. T* MSA = 0.07858 to 
be compared with the exact T* — 0.04917(2)2^, whereas the critical density is largely underestimated, i.e. p* c MSA = 
0.01449 to be compared with the exact p* = 0.080(5) 38 . By contrast the OMF scheme gives a slightly better critical 



0.12 




FIG. 3: Coexistence curves of the RPM in the one-loop approximation. From top to bottom: OMF scheme, WCA scheme, 
MSA scheme. Solid lines : coexistence curves, dashed line: spinodal curves. 

density than the MSA scheme - however still largely underestimated-, i.e. p* c OMF = 0.01834, but overestimates T* , i.e 
T* OMF = 0.1150. The results obtained in the framework of the WCA scheme lie somewhere between those obtained 
via the OMF and MSA schemes. 



2. Two-loop results 

The critical parameters of the RPM at the two-loop order are given in table JllJ and the coexistence curves are 
sketched in figure Including the two- loop correction yields a slight decrease of the critical temperature in all cases 
(except for the MSA case) as well as an increase of the critical density. More precisely one finds, at the two-loop level, 
the following critical parameters : (T* OMF = 0.1133, p* OMF = 0.0316), (T* MSA = 0.07993, p* MSA = 0.01722), 
and (T* WCA = 0.08428, p* WCA = 0.0137) in the OMF, MSA, and WCA schemes respectively. The improvement 
upon one-loop results is therefore significant but the theoretical predictions, notably those for the critical densities, 
are still far from the MC results. Our theoretical results can be (unfavorably) compared with those obtained recently 
in the framework of the the so-called collective variable (CV) method developed by the Ukrainian school which yield 
{T* = 0.0502, p* = 0.042)22i2i. At the moment the relation between the CV method and the KSSHE formalism is 
not established, work in this direction in under way. 

At the two-loop order the free energy of the binary SPM depends explicitely on the charge assymetry factor £. The 
effects of £ on the coexistence curve parameters in the OMF, WCA, and MSA schemes are resumed in table Ijlllll for 
£ = 1,2,3,4,5. In all the considered regularization schemes the critical temperature T* is found to increase slightly 
with £ in clear contradiction with MC simulation data where a rapid decrease of T* with £ has been observed^. 
However, as apparent in table Ijllfr. the critical densities -despite much too low values- are found to increase slightly 
with £ in the OMF and MSA schemes, which is in qualitative agreement with the MC data of ref. (|34j). 
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FIG. 4: Coexistence curves of the RPM in the two-loop approximation. From top to bottom: 
MSA scheme. Solid lines : coexistence curves, dashed line: spinodal curves. 



OMF scheme, WCA scheme, 



V. SUMMARY 



In this paper we have studied the liquid-vapor coexistence curve of the SPM -i.e. a simple version of the primitive 
model, the prototype of a system governed by long range Coulomb interactions-, by means of various mean-field 
theories. These theories were obtained in the framework of the KSSHE field-theoretical representation of charged HS 
systems by means of a loop expansion of the free energy of the homogeneous systeirJJ». Some of these MF theories 
are equivalent to well known approximations of the theory of liquids such that the Debye-Hiickel, or MSA (ORPA) 
theory, other were considered here for the first time. The results obtained for the critical point are in all cases in poor 
agreement with the MC simulations. 

The peculiarity of ionic criticality is clearly seen in the KSSHE formalism. The zero-loop order of the KSSHE 
theory is unable to reproduce a liquid- vapor transition and provides merely the local charge neutrality condition 11 ; 
the transition appears only at the one-loop order which is an unusual feature and suggests some kind of inconsistence 
of the KSSHE theory for ionic systems. The same remark holds for the other field-theoretical treatments of the 
SP M 19 ' 20 ' 21 ' 22 ' 24 . From the KSSHE side, further challenges are the derivation of the one-loop corrections to the 
density correlation functions and generalization of the two-loop expansions of the homogeneous free energy to the 
case of a PM with size assymetry. 
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APPENDIX A: THE RELATIONS BETWEEN g£° t AND T 



We shall denote by V e ] the GCPF of the SPM in the presence of an external electrostatic potential 14(1). 

The M dimensionless chemical potentials v a are assumed to be uniform. Let /o e (l) be the external charge distribution 
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giving rise to V^(l). V c (l) and p e (l) are linked by the Poisson equation which can be written, with the help of the 
formal notations of sec. {I1J, as 

«- 1 (l ) 2)F e (2)=p e (l). (Al) 
Obviously the charge correlation functions defined at equation (|2.17|) can be expressed as the functional derivatives 

(-1)" 1 S n E[{v a },V e ] 



G<g\l,...,n) 



/?« S[{u a }] SV c (l)...SV c (n) 



v Q =o 



r (n)T n , _ (-1)" ^ln5[K},y c 

Lr c (±,...,71) - TtTTiI H77T 



/?" ^ (l)...^ (n) 



(A2) 



Ve = 



Expressions similar to equations l|A2jl can be derived for the field correlation functions and G^ T in the 

following way. Firstly note that one has S[{^ a }, V e ] = S[{^ Q — q a V c }} where the "smeared" field V c is given by 
V e (l) = /3t(1, 2) V c (2). Therefore the expression (|2.13|l of S can be written more explicitely as 

S[{u a },V ] =KJ-J exp(-i(^| U - 1 |^)-lnS H s[{^ + z ( z ce (0-F c )}]). (A3) 

We now perform the change of variables 

ip -> if - ip^ 2 V e (A4) 
in equation (|A3|) . The functional Jacobian is equal to one and this gives us, after taking the logarithm 

haS[K},K] = ~ (V^lVe) +]nZ*[{v a },B] , (A5) 
where lnS*[{V Q }, B] is the KSSHE partition function l|2.13[) associated with the modified action 

(A6) 

where -B(l) = i(3 1 / 2 v~ 1 (l : 2)V c (2) plays the role of an additional external magnetic field. As it is well known, the 
functionals S* [-{V Q }, B] and ln2*[{^ , Q ,},i3] are the generators of the ordinary and connected correlation functions of 
the KSSHE field ip respective!}* 1 ^. Therefore we have 



G^(l, 
G< n)T (l, 



,«-) = - 



l 5 n E*[{u a },B] 



Z[{v a }] SB(l)...6B(n) 
5 n \nS*[{v a },B] 



B=0 



SB{l)...5B(n) 



B=0 



Taking the functional derivative of both sides of equation (|A5|) with respect to V c and noting that 
one finds that 



pc(l)+Pe(l) = -*P- 1/2 V- I (l,2)<<p(l)> 



(A7) 



(A8) 



(A9) 



~ 4^/2 

yielding the charge neutrality condition pc(l) — for a uniform system in the absence of external charge, since in 
that case < <p(l) >= cte . 

Taking twice the functional derivative of both sides of equation (|A5|) with respect to V e in the limit B — > yields 
readily 



(A10) 



PG {2) T (l, 2) = ^A 1( K1, 2) - JL_A! A 2 G^ T (l, 2) , 

where we made use of equations (|A2I . IJA7I) . and l|A8() . Finally, continuing the process and differentiating equation (|A5|) 
n times (n > 3) one obtains 



n/2 -n G (n) _ _ ^ = L-^ Al . . . A„G^ n) T (l, ...,«) 

Equations i|A9|l . (|A10|) . and (|A11|) were obtained in a complicated way in ref. |ll|) . 



(All) 
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APPENDIX B: AN ALGORITHM TO DETERMINE THE COEXISTENCE CURVE 



In this appendix we give some details on the algorithm used to obtain the coexistence curve of the SPM in the 
various approximation schemes considered in the paper. We denote by w^""' (v, (3, (p) the n-th partial derivative of 
the Landau function with respect to the order parameter ip. For n > 2, w ( n ) (v, /3, ip) is independent of the chemical 
potential v as apparent from equation (|4.1() . At any temperature T the function w^ 3 '(P,ip) versus <p has the shape 
depicted in figure i.e. with a unique zero located at some ip$. A precise numerical estimate of tp^ can be obtained 
with the help of the bisection method. In that aim we used the routine RTBIS of Numerical Recipes 40 -. More precisely 

<p 3 = RTBIS( W ( 3 \^ min ,(^ max ,e) , (Bl) 

where RTBIS seeks for the (unique) zero of w^ 3 ' (ip) in the prescribed interval [</? m i n , <£> m ax] with some arbitrary precision 
e. 

Therefore the function w^'{j3, ip) versus <p is convex for all (3. For T < T c has the shape depicted in figure QJ, 
i.e. with a unique negative minimum at tp — ip^. For T > T c the minimum is positive and thus the critical temperature 
T c is obtained by the condition (/3 C , ipz) = 0. The zeros of w^(ip) for T < T c will be denoted ipg and ipf > cpg. 
They are also easily obtained with the help of the bisection method, i.e. 

ip$ = KTBlS(w^\ip min ,ip 3 ,e) , 

<pf = RTBIS( W ( 2 \^ 3 ,Vmax,e) • 

Since there the second derivative of the free energy with respect to p (i.e. the density) vanishes, the locus of the 
points (iPgipf) is the so-called spinodal curved. 

The points (ff,ff) correspond to the maximum and the minimum of the function w^'(<p) respectively. For an 
arbitrary chemical potential v 7^ v C ocx one is likely to find the dashed curves of (ip) and w{ip) displayed in figure 
One wants to find the coexistence chemical potential is CO cx such that the two minima of w(ip) have the same value, i.e. 
to find for w((p) and w^(p) the solid curves depicted in figure (JIJ. This can be achieved conveniently and precisely 
by means of the following iterative process. 

(i) First one needs an estimate of v not too far from v cac ^. Denoting by «J~ in and Wmlx the values of the minimum 
and maximum of w^(ip) respectively one adopts for v the new value v — > vq = v + (w^. + WmL)/2 which ensures 

+1, j- (!) (!) 

that nOW W^( n = -Wmax- 

(ii) The zeros of w^(ip) in the intervals [<£> m in, Pg] & n d [pf , V?max] are rough estimates of the coexistence densities. 
Denoting them by p° g and pf respectively one has 

<p° g = RTBIS^yw^.e) , 
= KTBlS(w < - 1 \ pf , ip max , e) . 

(hi) At this point we have found a z/o and two estimates pP g and p® such that u>W ( l / 0) ip°) = w ( x ) (;y 0; y>°) = but 
w(vo, (p° g ) 7^ w{vq, ip®) in general. It follows from equation 1)4. l|l that if one defines 

w{vo,f°i) - w(is ,p°) 
+ 5 — , B2 

<Pl - <Pg 

then w(vi, ip° g ) — w(ui, ip®) but w^'^yi, p° g ; ) 7^ 0. One then determines two new estimates of ipi and (p g as 

p\ = RTBIS(w (1) (^i,< y 5),^ mi n,V5f ,e) , 
<p\ = KrBIS(w^(vi,ip),tpf ,tp max ,e) . 

Now one has w^'(fx, ipg) — w^ l \ui, p\) — and one computes again w(y\, ip g ) and w(i>\, pj). If 



Mvi,<p g ) -w{vi_,(p}) 



< e , (B3) 



where e is the wanted precision, then the problem is numerically solved. If it is not the case then one sets vi — > vq 1 
LPg~ j — > ipPg i and one goes back to point (iii). For any reasonable e a few iterations are necessary to get the result with 
a precision e. 
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We end this appendix by a few comments on the numerical calculation of w and its derivatives w^ n \ The HS free 
energy and its derivatives were computed in the Carnahan-Starling approximation^. We also tried the PY expressions 
(in the compressibility route) which makes little difference at the rather low densities considered here. Only the 
numerical estimate of the functions K[X, p] (cf equation l|3.5af> ) which enter the free energy at the second-loop order 
needs additional comments. In order to compute K[\, (p] and its partial derivatives with respect to ip (note that 
A = 47r flip does also depends on tp) we have used the PY solution for has, tp- For numerical purposes it is necessary 
to reexpress K in terms of the function 7hs, <p{r) = has, tp( r ) ~ chs. u>(?) which is a continuous function of r. In the 
framework of PY theory analytical expressions can be obtained for the Fourier transform 7hs, tp and its derivatives 
with respect to tp with the help of a computer algebra package. 7hs, tp( r ) & n d its derivatives with respect to <p are 
then obtained as inverse numerical Fourier transforms. 
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TABLE I: One-loop results for the critical temperature T c , chemical potential v c , pressure P c , density p c , and energy u c of the 
SPM. Various regularization schemes are considered, as explained in the text. 





T c 


v c 


10 2 x I3P C 




10 x u c 


DH 


0.3271 


-5.5274 


2.4194 


0.12267 


-1.3316 


T a : a* = 0.1 


0.2868 


-5.8553 


1.8336 


0.10015 


-0.9552 


a* = 0.2 


0.2532 


-6.1688 


1.4052 


0.08224 


-0.6981 


a* = 0.3 


0.2250 


-6.4703 


1.0860 


0.06778 


-0.5170 


a* = 0.4 


0.2012 


-6.7608 


0.8453 


0.05598 


-0.3866 


a* = 0.5 


0.1809 


-7.0408 


0.6622 


0.04630 


-0.2913 


a* = 0.6 


0.1636 


-7.3104 


0.5220 


0.03835 


-0.2210 


a* = 0.7 


0.1487 


-7.5695 


0.4141 


0.03180 


-0.1686 


a* = 0.8 


0.1358 


-7.8182 


0.3307 


0.02641 


-0.1293 


a* = 0.9 


0.1247 


-8.0564 


0.2659 


0.02198 


-0.09977 


OMF 


0.1150 


-8.2843 


0.2154 


0.01834 


-0.07745 


WCA 


0.08446 


-9.1737 


0.0930 


0.00880 


-0.02813 


MSA 


0.07858 


-9.1393 


0.1229 


0.01449 


-0.04864 



TABLE III: Two-loop results for the critical temperature T c , chemical potential v c , pressure P c , density p c , and energy u c of 
the assymetric SPM. Various regularization schemes are considered, as explained in the text. 





T c 




10 2 x f3P c 




10 X Uc 


OMF £ = 1 


0.1133 


-7.8229 


0.4181 


0.0316 


-0.1363 


OMF £ = 2 


0.1263 


-7.0281 


0.5848 


0.0336 


-0.1209 


OMF £ = 3 


0.1433 


-6.1618 


0.8049 


0.0367 


-0.0966 


OMF £ = 4 


0.1585 


-5.5175 


1.0011 


0.0396 


-0.0691 


OMF £ = 5 


0.1719 


-5.0252 


1.1744 


0.0421 


-0.0399 


WCA £ = 1 


0.08428 


-8.9142 


0.1431 


0.0137 


-0.0461 


WCA £ = 2 


0.09334 


-8.2044 


0.1903 


0.0132 


-0.0375 


WCA £ = 3 


0.1058 


-7.4035 


0.2582 


0.0137 


-0.0305 


WCA £ = 4 


0.1171 


-6.8020 


0.3208 


0.0146 


-0.0244 


WCA £ = 5 


0.1272 


-6.3422 


0.3770 


0.0155 


-0.0183 


MSA £ = 1 


0.07993 


-9.0623 


0.1379 


0.0172 


-0.0602 


MSA £ = 2 


0.08422 


-8.4828 


0.1926 


0.0199 


-0.0628 


MSA £ = 3 


0.09060 


-7.7191 


0.2751 


0.0231 


-0.0602 


MSA £ = 4 


0.09677 


-7.0701 


0.3555 


0.0257 


-0.0527 


MSA £ = 5 


0.1025 


-6.5333 


0.4302 


0.0279 


-0.0422 



TABLE II: Two-loop results for the critical temperature T c , chemical potential v c , pressure P c , density p c , and energy u c of 
the RPM. Various regularization schemes are considered, as explained in the text. 





Tc 


Vc 


10 2 x I3P C 


Pc 


10 x u c 


OMF 


0.1133 


-7.8229 


0.4181 


0.0316 


-0.1363 


WCA 


0.08428 


-8.9142 


0.1431 


0.0137 


-0.0461 


MSA 


0.07993 


-9.0623 


0.1379 


0.01722 


-0.06021 



